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Abstract. We present a comparison of different numerical techniques for the 
integration of variational equations. The methods presented can be applied 
to any autonomous Hamiltonian system whose kinetic energy is quadratic in 
the generalized momenta, and whose potential is a function of the generalized 
positions. We apply the various techniques to the well-known Hcnon-Hcilcs 
system, and use the Smaller Alignment Index (SALI) method of chaos detection 
to evaluate the percentage of its chaotic orbits. The accuracy and the speed 
of the integration schemes in evaluating this percentage are used to investigate 
the numerical efficiency of the various techniques. 

1. Introduction. The determination of the stability of motion is of great impor- 
tance when investigating nonlinear dynamical systems. To distinguish correctly 
between regular and chaotic motion, several different methods have been developed 
during the years. Most of these techniques, such as the maximal Lyapunov expo- 
nent [15], the fast Lyapunov indicator [4] or the Smaller Alignment Index (SALI) 
[14] , rely on the study of the time evolution of deviation vectors from a given orbit 
to discriminate between the two regimes. The time evolution of these vectors is 
governed by the so-called variational equations. 

Besides the correct determination of the regular or chaotic nature of individual 
orbits, in many cases, statistical statements over a large region of the phase space 
are also needed. For example, in order to determine the percentage of regular and 
chaotic orbits in a given system, the characterization of many orbits is required. 
In addition to the accurate computation of chaos indicators also for such tasks the 
CPU time needed to perform these computations becomes very important. In this 
work we compare different numerical techniques for the integration of the variational 
equations, concentrating on their accuracy and computational speed. 

The paper is organized as follows: in section 2 we describe the general layout of 
our investigation. We concentrate our study on a simple two degrees of freedom 
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Hamiltonian system: the well-known Henon-Heiles system [7], which is presented 
in section 2.1. In our study we use the SALI method of chaos detection, which 
is presented in section 2.2. To solve the equations of motion of the Hcnon-Hcilcs 
system and the associated variational equations one has to employ numerical meth- 
ods. Any non-symplectic general-purpose integrator can be used for this task. In 
sections 2.3.1 and 2.3.2 we present two such techniques, which we use in our study. 
In [16] it was shown that it is also possible to use methods based on symplectic 
integration techniques to solve these equations. In section 2.3.3 we describe shortly 
the most efficient of these techniques: the so-called tangent map (TM) method. A 
numerical procedure to obtain relatively fast information on the nature of orbits for 
a large set of initial conditions is then given in section 2.4. In section 3 we present 
our numerical results for individual orbits, as well as global results for the whole 
system. The summary and the conclusions of our study are found in section 4. 

2. Numerical integration of variational equations. 

2.1. Henon-Heiles system. The Hamiltonian function of the Henon-Heiles sys- 
tem [7] is 

H(q 1 ,q 2 ,p 1 ,p 2 ) = \{p\+pD + ^(<7? + <ll) + 9??2 - ^f, ( l ) 

with qi, q 2 being the generalized coordinates and pi, pi the conjugate momenta. 
The orbit evolution is given by the standard Hamilton equations of motion 

9H , . dH 

qt = ^— and pi = - — , « = 1,2, (2) 
dpi oqi 

where the dot denotes derivation with respect to time t. The time evolution of the 
variations Sqi , Spt (which can be considered as coordinates of a deviation vector) is 
governed by the variational equations, given by 

2 d 2 H 

Sqi = Spi and 5p q = - ^ <%, i = l,2. (3) 

j=i vqiCqj 

Eqs. (2) and (3) form a coupled system of ordinary differential equations. It should 
be noted that the solution of (3) depends explicitly on the solution of (2), i. e. on 
the reference orbit qi(t),Pi(t), and thus Eq. (3) cannot be solved independently from 
Eq. (2). 

2.2. SALI method. The evaluation of the SALI is an efficient and simple method 
to determine the regular or chaotic nature of orbits in dynamical systems. The 
SALI was proposed in [14] has since been successfully applied in order to distinguish 
between regular and chaotic motion both in symplectic maps and Hamiltonian flows 
[17, 18, 13, 12, 3, 11, 19]. For the computation of the SALI of a given orbit, one 
has to follow the time evolution of the orbit itself and also of two deviation vectors 
Vi(t), V2(t), which initially point in two different directions. Then, according to [14] 
the SALI is defined as 

SALI(t) = min + V 2 (t)j, ||ft(t) - ? 2 (t)||} , (4) 

where || • || denotes the usual Euclidean norm and Vi, i = 1,2 are normalized vectors 
with norm equal to 1. 

The SALI has a completely different behavior for regular and chaotic orbits, 
and this allows us to clearly distinguish between them. In particular, the SALI 
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fluctuates around a non-zero value for regular orbits, while it tends exponentially 
to zero for chaotic orbits [14, 17], following a rate which depends on the difference 
of the two largest Lyapunov exponents [18]. 

2.3. Used numerical methods. 

2.3.1. DOP853. The DOP853 1 integration method belongs to the big class of ex- 
plicit Runge-Kutta methods. This non-symplectic scheme of order 8 is based on the 
method of Dormand and Price (see [5, Sect. II.5]). We use this integrator to solve 
the set of differential equations composed of Eqs. (2) and (3). Two free parameters, 
r and 8, are used to control its numerical performance. The first one defines the 
time span between two successive outputs of the computed solution. After each step 
of length r the deviation vectors are renormalizcd and the value of SALI is com- 
puted. For the duration of each step r. the integrator adjusts its own internal time 
step in order to keep the local one-step error smaller than a user-defined threshold 
value 6. For the DOP853 integrator the estimation of this local error and the step 
size control is based on embedded formulas of orders 5 and 3. 

2.3.2. Taylor methods. The basic idea of the so-called Taylor series methods (for 
details see for example [5, Sect. 1.8] and references therein) is to approximate the 
solution at time U + t of a given s-dimensional initial value problem 

^ = /(!/(*)) yeR s ,teM (5) 

dt 

from the nth degree Taylor series of y(t) at t = ti as 

, . dy(U) T 2 d 2 y(tA t" d n y(t t ) 

The computation of the derivatives is commonly done using automatic differentia- 
tion (see for example [8]). 

In our study we use two different public available implementations of the Taylor 
method: TIDES 2 [2] and TAYLOR 3 [8]. Both methods have internal automatic 
order and step size computation to ensure the user-defined local one-step error S. 
Also here the parameter t defines the step size, after which the renormalization of 
the deviation vectors and the computation of SALI is done. 

The whole testbed of our work is written using the FORTRAN programming 
language exploiting extended double precision 4 . While TIDES offers directly a 
FORTRAN integration routine, a wrapper to include the routine written in C had 
to be used for TAYLOR. Therefore for the latter only 16 significant digits were 
available for the integration. 

2.3.3. TM method. Besides general-purpose integrators, it is also possible to use 
techniques based on symplectic methods to integrate the Hamilton equations of 
motion and the corresponding variational equations. This was shown in [16], where 
a thorough discussion of possible methods can be found. The most effective of 
these techniques, the TM method, is used in this work. Let us outline its basic 
idea, which is founded on a general result stated for example in [9]: Symplectic 
integrators can be applied to first order differential systems X = LX, that can be 



1 Freely available from http://www.unige.ch/~hairer/software.html. 
2 Freely available from http://gme.unizar.es/software/tides. 

3 Freely available from http://www.maia.ub.es/~angel/taylor/software/. 

4 Corresponding to 18 significant digits or equivalently to a machine accuracy of ~s 10" 
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written in the form X = (La + Lb)X, where L,La,Lb are differential operators 
defined as L x f = {%, /} and for which the two systems X = LaX and X = LbX 
are integrable. Here {/, <?} are Poisson brackets of functions f(q, p), g(q, p) defined 
as: 

The set of Eqs. (2) and (3) is one example of such a system, because Hamiltonian 
(1) can be divided into two integrable parts A and B with H = A(p) + B(q). A 
symplectic integrator splits the equations of motion (2) into several parts, applying 
either the operator La or Lb ■ These are the equations of motion of the Hamiltonians 
A and B, which can be solved analytically, giving explicit mappings over the time 
step CjT, where the constants c, are chosen to optimize the accuracy of the integrator. 
These mappings can then be combined to approximate the solution after time step 
r. In [16] it was shown that the derivative of these mappings - with respect to 
the coordinates and momenta of the system (the so-called tangent maps) - can 
be used for the time evolution of deviation vectors or, in other words, for solving 
the variational equations (3). We note that the TM method is called the 'global 
symplectic integrator method' in [10]. 

In [9] a family of symplectic integrators called SABA„ and SBAB n was intro- 
duced, with n being the number of applications of operators La and Lb- These 
integrators have only positive intermediate steps and can be used with an addi- 
tional corrector step C at the beginning and the end of each step r to increase their 
accuracy. An integrator of order 4 of this family, namely the SBAB2C integrator 
which includes corrector steps, is used in our investigation. A detailed description of 
the application of the SBAB2C integrator for the TM method to the Hcnon-Hcilcs 
system can be found in [16]. 

2.4. Fast PSS method. Besides information on the chaotic or regular character of 
individual orbits, a more global description of dynamical systems is also of interest. 
For example, such a study could include the computation of the percentage of 
regular/chaotic orbits for a given set of initial conditions (ICs). This information 
requires the integration of the equations of motion/ variational equations, and the 
computation of a chaos indicator for the whole set of ICs, which can become a very 
hard computational task. In order to address this problem, we implement a method 
proposed in [1], which exploits the Poincare surface of section (PSS) of the system 
in order to speed up this computation. 

For a fixed value of Hamiltonian (1) (throughout our study we use always H = 
0.125) we define the PSS of the system as the plane given by q\ = and pi > 0. 
Each point in this plane defines a set of values (92,^2)- To evaluate the percentage 
of regular orbits, one normally computes for each point the value of some chaos 
indicator using a dense set of points on the PSS as ICs. In Fig. 1 we consider a grid 
of 400 x 400 ICs and color each one according to its SALI value at t = 3000. 

Each orbit starting from any IC intersects the PSS in many points, and so, 
its SALI value can be attributed to all orbits having these intersection points as 
ICs. Therefore all these ICs do not have to be integrated separately. This procedure 
decreases drastically the CPU time needed for the global description of the system's 
chaoticity. We refer to this approach as the fast PSS method. 



3. Numerical results. 
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Figure 1. The PSS of the Henon-Heiles system (1) for qi = and 
pi > 0. A grid of 400 x 400 initial conditions is used and each 
initial condition is colored in grey-scale according to its SALI value 
at t = 3000. Black dots forming five closed curves correspond to 
the regular orbit Rl of Fig. 2, while the scattered black dots are 
intersection points of the chaotic orbit CI of Fig. 3 with the PSS. 



3.1. Individual orbits. Let us first investigate, how well the different methods 
described in section 2.3 can determine the nature of individual orbits. We use these 
methods to integrate Eqs. (2) and (3), and then we compute the evolution of the 
SALI in order to determine the nature of the orbit. Unless otherwise stated, we 
always renormalize the deviation vectors after each time step of length r = 0.05. 
For the non-symplectic routines we adopt a one-step accuracy of S = 10~ 5 . 

As representative examples we consider 3 orbits of the Henon-Heiles system with 
different dynamical behaviors. The evolution of the SALI for a regular orbit (Rl) is 
presented in Fig. 2, while in Fig. 3 we have similar results for a chaotic orbit (CI). 
Finally, in Fig. 4 the SALI of a sticky chaotic orbit (C2) is shown. 




Figure 2. Time evolution of the SALI for the regular orbit Rl 
with initial conditions q\ = 0, pi « 0.2334, q2 = 0.558, P2 = 0. 



From Fig. 2 we see that the results obtained for orbit Rl by the different in- 
tegration methods are nearly identical. As theory predicts a constant SALI for 
regular motion, such behavior is correctly identified for orbit Rl. Information con- 
cerning the numerical performance of various techniques for orbit Rl is given in 
Table 1 (throughout this paper the reported CPU times refer to an Intel Xeon 
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Figure 3. Time evolution of the SALI for the chaotic orbit CI 
with initial conditions q\ = 0, p\ « 0.4208, q% = —0.25, P2 = 0. 
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Figure 4. Time evolution of the SALI for the sticky chaotic orbit 
C2 with initial conditions q\ = 0, p\ « 0.11879, (72 = 0.335036 and 
p 2 = -0.385631. 



X5660, 2.80 GHz computer). From the results reported in this table, we see that 
differences between the applied methods appear in their energy conservation prop- 
erties, as is indicated by the relative energy error \AH/H\, shown in Fig. 5. For 
the symplectic algorithm (TM method) \AH/H\ shows fluctuations around 10~ 8 , 
while it grows with time for the other methods as expected (see for example [6, 
Sect. IX.8]). The TAYLOR method has the worst performance since it is able to 
conserve the energy only up to an error level of ~ 10~ 6 for S = 10" 5 . This is 
probably due to the 2 digits less in accuracy that are available for this method 
(see section 2.3). The best method with respect to the energy conservation is the 
TIDES algorithm for which the relative error is w 10~ 13 (S = 10~ 5 ) and « 10~ 16 
(5 = 10~ 16 ) at t = 10 7 . The price paid for the excellent accuracy of the algorithm 
is that TIDES requires, in general, the largest CPU times and the highest orders 
among the tested methods. 

Orbit CI is also correctly identified as chaotic by all methods in less than 1000 
time units, within which SALI goes to zero (Fig. 3). A difference is found in the 
results for the sticky chaotic orbit C2 (Fig. 4). Up to t w 10 6 all methods indicate 
a regular behavior. It is only afterwards that SALI goes to zero for the TM, the 
DOP853 and the TIDES methods, correctly identifying C2 as a sticky chaotic orbit. 
For the used values of S and r the TAYLOR method does not succeed to show the 
decrease of SALI to zero, at least up to t = 10 7 . 

3.2. Global results. In order to find the percentage of regular and chaotic orbits 
of the Henon-Heiles system (1), we compute for a grid of 400 x 400 ICs on the PSS 
the SALI value for different final times tfi na i. One could argue that due to the finite 
resolution with which the grid of ICs is taken on the PSS, the fast PSS method 
(Sec. 2.4) would not be as accurate as the individual computation of SALI for each 
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Table 1 . Information on the performance of the different numeri- 
cal methods used for the computation of the evolution of the regular 
orbit Rl, of two deviation vectors from it and of its SALI. The step 
size r was always 0.05. The order used by the TIDES and TAY- 
LOR methods is determined by these routines for each step r and 
is constant for the whole integration. 
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Figure 5. The time evolution of the relative energy error \AH/H\ 
for the different integration schemes in the case of the regular orbit 
Rl. The step size r was 0.05 for all methods. For all non-symplectic 
methods we set (a) S = 10 -5 and (b) 5 = 10~ 16 . The relative 
energy error for the TM method is the same for both panels and 
is reported only for reference, since for this method no one-step 
accuracy S is defined. 

IC. The percentages (over the total number of ICs compatible with Hamiltonian 
(1)) of regular (SALI > 10~ 4 ), chaotic (10~ 8 > SALI), and sticky chaotic orbits 
(10~ 4 > SALI > 10~ 8 ) obtained by both approaches using the TM method, are 
shown in Fig. 6(a), while the required CPU times are reported in Fig. 6(b). From 
the results of Fig. 6 we see that both approaches obtain practically the same values, 
while the CPU time needed by the fast PSS method remains considerably smaller 
with respect to the full integration of individual orbits. For this reason we apply 
the fast PSS method for computing the percentages of regular, chaotic and sticky 
chaotic orbits for different values of the time step r and the final time ifinai- The 
obtained results can be found in Table 2. 

From the results of Table 2 we see that for large values of r and ifmai (t = 0.50 
and tfinai = 10 6 ) the non-symplectic methods find 3 — 4% less regular orbits than 
the TM method. In order to understand this discrepancy we computed for orbit Rl 
the evolution of the SALI by the DOP853, the TIDES and the TAYLOR methods 
with r = 0.50 and 5 = 10~ 5 (Fig. 7). From Fig. 7 we see that all non-symplectic 
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Figure 6. (a) Percentage of regular, chaotic and sticky chaotic 
orbits as a function of tfi na i, when each initial condition is inte- 
grated until time tfi na i (black dots) and when the fast PSS method 
was used (solid lines). The orbits are characterized according to 
their SALI value at time t& na \ as regular (SALI > 10~ 4 ), chaotic 
(1(T 8 > SALI) and sticky chaotic (1(T 4 > SALI > 1(T 8 ). A 
grid of 400 x 400 initial conditions on the PSS of Fig. 1 was used. 
The integrations of the orbit and the deviation vectors were done 
by the TM method using the SBAB2C integrator with a step size 
r = 0.05. (b) The CPU time needed for the computation of the 
results shown in (a). 



methods fail to detect correctly the regular character of the orbit because SALI 
drops to zero after t = 10 6 . Such behaviors lead to the increase of the percentage of 
sticky chaotic orbits in Table 2, since some regular orbits are wrongly characterized 
as sticky or chaotic. Decreasing S to values < 10 -14 solves the problem, as it leads 
to a correct identification of the orbit's nature, but also increases the required CPU 
time (see Fig. 7). 
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Figure 7. Time evolution of the SALI for the regular orbit Rl 
for r = 0.50. For the DOP853, the TIDES and the TAYLOR non- 
symplectic methods we set S = 10~ 5 (black line) and 5 = 10~ 14 
(grey line). The required CPU times are given as labels of the 
various lines. 



For smaller values of r the results obtained from different techniques are more 
consistent. For large integration times all methods give a similar percentage of 
regular orbits of 40%. TM, DOP853 and TIDES agree here within 0.01%. For 
t = 0.01 and t = 10 6 the TM method clearly discriminates between regular and 
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Table 2. Percentage of regular, chaotic and sticky chaotic orbits 
for different values of the final time ta na \ and the step length r, for 
the integration schemes presented in section 2.3. A grid of 400 x 400 
initial conditions on the PSS of Fig. 1 is used. 8 was set to 10 -5 for 
the non-symplectic methods. The required CPU times are reported 
in the last column. 
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method 


% regular 


To sticky chaotic 


7o chaotic 


CFU time 
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TM 
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h 02 min 






DOP853 
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h 04 min 
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0.50 
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0.07 
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h 09 min 
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0.37 
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h 17 min 
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h 15 min 
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0.00 
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DOP853 


36.94 


2.59 


60.47 


1 h 29 min 






TIDES 


35.78 


3.71 


60.51 


2 h 06 min 






TAYLOR 


34.01 


6.23 


59.76 


1 h 15 min 


0.01 


10 4 


TM 


40.38 


0.72 


58.90 


h 19 min 






DOP853 


40.22 
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h 38 min 
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0.84 
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TM 


39.39 


0.57 


60.04 


2 h 02 min 
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0.22 
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chaotic orbits, while all other methods still find w 0.3% of sticky chaotic orbits. But 
since the SALI of those orbits will eventually go to zero as well, it can be expected 
that also the non-symplectic methods will give the same results as the TM method 
when the integration time is increased. 

4. Summary and conclusions. We considered the problem of fast and accurate 
integration of the variational equations of a conservative Hamiltonian system. We 
compared different numerical techniques for this task, and applied them to the 
Henon-Hciles system. We considered non-symplectic methods of high accuracy; 
particularly the DOP853 scheme, as well as the TIDES and TAYLOR packages, 
which arc based on Taylor expansion techniques. We also applied the TM method, 
which exploits symplectic integrators, using in particular the SBAB2C integrator. 

The variational equations govern the evolution of small deviations from a given 
orbit. Using the SALI chaos indicator, which is defined through the time evolution 
of deviation vectors, we determined the chaotic or regular nature of individual 
orbits. In addition, applying an efficient numerical approach, the so-called fast PSS 
method, we were able to rapidly identify regions of order and chaos in the phase 
space of the system. 
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Our numerical results show that the TM method had the best numerical perfor- 
mance both in accuracy and in speed, especially for large integration steps, when 
the other non-symplcctic schemes failed to compute accurately the fraction of reg- 
ular and chaotic motion. For moderate integration steps all applied methods gave 
practically the same results, with the TM method being always faster. 

Among the non-symplectic algorithms the TIDES was the most accurate one 
producing similar results as the DOP853 integrator. In many cases the results of 
the TAYLOR method were found to be less accurate than the ones obtained by the 
other methods, probably due to some implementation peculiarities of the algorithm. 
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